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Abstract. We study the 6- vertex model with fixed boundary conditions. In the thermodynamical 
limit there is a formation of the limit shape. We collect most of the known results about the 
' analytical properties of the free energy of the model as the function of electric fields and study the 

^SJ . asymptotical behavior near singularities. We also study the asymptotic of limit shapes and the 

structure of correlation functions in the bulk. 
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1. Introduction 

Ising and dimer models were among the first models in two-dimensional statistical mechanics 
for where the partition function and for some of the correlation functions were computed explicitly 
in terms of Pfaffians of certain matrices. For this reason both of these models can be regarded as 
theories of Gaussian discrete two dimensional fermi field. The Ising model was solved by Onsager 
and for dimer models the Pfaffian solution was found by Kasteleyn. 

The 6-vertex model, a particular case of which is the ice model is interesting for a number of 
reasons. Physically, it is a model of ferro- and antiferro- electricity. It has many equivalent reformu- 
lations, one of them (which we will use) describe the 6-vertex configurations on a planar connected 
simply connected region in terms of stepped surfaces. One of the combinatorial reformulations of 
the 6-vertex model for specific value of parameters |Kup| is related to alternating sign matrices 
|Bres| . 

The 6-vertex model generalizes dimer models and can be regraded as the theory of Gaussian 
discrete fermions with four fermionic interaction. The partition function of the 6-vertex model 
with periodic boundary conditions was computed in [Lieb ] using the Bethe Ansatz method. 

The computation of correlation functions in the 6-vertex model is highly non-trivial and to the 
large degree is still a challenge. There are two known approaches to this problem based on the 
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internal symmetry of the model, i.e. on the representation theory of quantum affine algebras. First 
approach is based on determinantal formulae for certain matrix elements |KBI| . For some recent 
results based on this method see |BPZ| and [CP] , The second approach is based on form-factor 
formulae derived in [Smj . For an overview of this approach see |JM| and for the latest results see 
|BJMST| . 

The 6-vertex model on a planar simply regions can be reformulated as the theory of random 
stepped surfaces. A configuration of arrows in a 6-vertex model can be interpreted as a configuration 
of paths which can be viewed as level corves of a height function defining a stepped surface. Gibbs 
measure on 6-vertex configurations define a Gibbs measure on stepped surfaces. 

For certain class of such Gibbs measures, random surfaces in the thermodynamical limit develop 
the limit shape phenomenon [Shef] also known as the acric circle phenomenon [CEP]. It means 
that on macroscopical scale the random surface becomes deterministic. The fluctuations remain 
at smaller scale and the structure of fluctuations may change depending on how singular the limit 
shape is at this point. The limit shape phenomenon is studied in details in dimer models |K0) . 

The numerical results from |AR]|SZ| show how limit shapes develop in a 6-vertex model with 
domain wall boundary conditions. In this paper we will focus on the limit shape phenomenon for 
the 6-vertex model on a planar connected simply connected regions with fixed boundary conditions. 
Computing the limit shape involves two steps. 

Step one is the derivation of the formula for the free energy of the 6-vertex model as a function of 
magnetic fields. Unlike the partition function of dimer models, the free energy can not be computed 
explicitly. However, it can be written in terms of the solution to a linear integral equation. Many 
analytical properties of the partition function are known but scattered in the literature, see for 
example |Yang| [YangT] [SY] [SYl] [W] [BS] [Nold] [NK] . We collected most of them in the section [3] 
together with some new results. 

Step two is to derive and solve the variational principle which determine the limit shape for 
given boundary conditions. The free energy of the model as a function of magnetic fields determine 
the functional in the variational problem. Such variational problem was first introduced for dimer 
models in [CKP| . The idea of using this variational problem in the 6-vertex first appeared in [ZlJ 
where some interesting partial results were obtained for correlation functions in the bulk of the 
limit shape. 

The structure of fluctuations near the limit shape is determined by the asymptotical behavior of 
correlation functions at smaller scales. We will discuss this problem for the 6-vertex model in the 
last section. 

Finally let us mention the special case of domain wall boundary conditions. These boundary 
conditions first appeared in the computation of norms of Bethe vectors [Korj . The remarkable fact 
about them is that the partition function of the 6-vertex model with these boundary conditions can 
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be written as a determinant |Izej . Another remarkable fact is that exactly these boundary conditions 
relate the 6- vertex model with alternating sign matrices |Kup| . The large volume asymptotic of the 
partition function of the 6- vertex mode with these boundary conditions was computed in [KZ][Z] . 

Here is the outline of the paper. In the first two sections we recall some basic facts about the 
6-vertex model, about its reformulation it in terms of height functions, and about the thermody- 
namical limit in the model. The third section contains the description of the free energy per site for 
as the function of electric fields in thermodynamical limit for periodic boundary conditions. Some 
asymptotical behaviors of the free energy are computed in section 4. This section is a combination 
of an overview and original results. In section 5 we study the asymptotical behavior of the limit 
shapes near the "freezing point". In the last section we discuss fluctuations. 

We thank C. Evans, R. Kenyon, A. Okounkov, and S. Sheffield for interesting discussions. The 
work of N.R. was supported by the NSF grant DMS 0307599, by the Niels Bohr initiative at Aarhus 
University, by the Humboldt foundation and by the CRDF grant RUMIl-2622 The work of K.P. 
was supported by the NSF RTG grant and DMS 0307599. 

2. The 6- Vertex Model 

2.1. The 6-vertex model. First, let us fix the notation. A square N x M grid L^^m is a graph 
with 4- and 1-valent vertices embedded into such that 4-valent vertices are located at points 
(n, m), n = 0, 1, . . . , — 1, m = 0, 1, . . . , M — 1 (see Fig. [2]) and 1-valent vertices are located 
at (-1, m), (A^, m), m = 0, 1, . . . , M - 1 and at (n, -1), (n, M - 1), n = 0, 1, . . . , A^ - 1. An edge 
connecting two 4-valent vertices is called an inner edge and an edge connecting a 4-valent vertex 
with a 1-valent vertex is called an outer edge. 

States of the 6-vertex model on L7v,m are configurations of arrows assigned to each edge (i.e. 
orientations of L^v,*/)- They satisfy the ice rule: at any vertex the number of incoming arrows 
should be equal to the number of outgoing arrows. Six possible configurations at a vertex are 
shown on Fig. [TJ Configurations of arrows on boundary edges are called boundary conditions. 

Each configuration of arrows on the lattice can be equivalently described as the configuration of 
"thin" and "thick" edges or empty and occupied edges shown on Fig. [H There should be an even 
number of thick edges at each vertex as a consequence of the ice rule. The thick edges form paths. 
We assume that paths do not intersect (two paths may meet at an ai-vertex). So, equivalently, 
configurations of the 6-vertex model can be regarded as configurations of paths satisfying the rules 
from Fig. [TJ 

To each configuration of arrows on edges adjacent to a vertex we assign a Boltzmann weight, 
which we denote by the same letters. The physical meaning of a Boltzmann weight is exp(— ^), 
where E is the energy of a state and T is the temperature (in the appropriate units). Thus, all 
numbers ai, 02, 61, ci, and C2 should be positive. 
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Choosing the scale such that T = 1, it is natural to write Boltzmann weights in the exponential 
form. 

Cl=e-^^ C2 = e-^^ 

where Ei, E2, and E3 are dimensionless interaction energies of arrows at different types of vertices, 
and H and V are dimensionless horizontal and vertical components of the magnetic field, respec- 
tively. In this interpretation arrows are spins interacting with the magnetic field. We set ci = C2 
because for the types of boundary conditions we will consider the difference between the number of 
ci-vertices and C2 vertices is the same for all states and, therefore, the probability does not depend 
on the ratio ci/c2. 

We also use the standard notation 

These are the weights of the model when there is no magnetic field. 

The weight of a state is the product of weights of vertices in the state. The weight of a state on 
L]\j (up to a constant factor) can be written in terms of energies and magnetic fields as 

exp(-^iAf(a) - E2N{b) - E^N{c) + ^N{hor) + ^A^(t;er)) 

where N[a) is the total number of a-vertices, N[b) is the total number of 6-vertices, N{c) is the 
total number of c-vertices, N{hor) is the total number of horizontal edges occupied by paths, and 
N{vert) is the total number of vertical edges occupied by paths. 

The partition function is the sum of weights of all states of the model 

Z = If (vertex), 

states vertices 

where ^(vertex) is one of the weights from Fig. [TJ 

Weights define the probabilistic measure on the set of states of the 6-vertex model. The prob- 
ability of a state is given by the ratio of the weight of the state to the partition function of the 
model 



P{state) 



n vertices ^ (vertex) 



Z 

This is the Gibbs measure of the 6-vertex model. 

Let us define the characteristic function of an edge e as 



aeistate) 



1, if e is occupied by a path; 
0, otherwise 
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Figure 1 . The 6 types of vertices and the corresponding thin and thick edges configurations 

A local correlation function is the expectation value of the product of such characteristic func- 
tions: 



((Tei(Te2..0-e„) = ^ P(state) (Te, (state) 
states i=l 

2.2. Boundary Conditions. 

2.2.1. Let us fix arrows (or, equivalently, think edges) on outer edges of L^^m- States in the 6- 
vertex with the same configurations of arrows on the boundary are called states with fixed boundary 
conditions. The difference between two such states can occur only at inner edges. 

The space of states with fixed boundary conditions is empty unless the boundary values satisfy 
the ice rule: the total number of incoming arrows on the boundary edges should be equal to the 
total number of outgoing arrows. In the path formulation this means that the number of paths 
through North and West boundaries should be equal to the number paths through the South and 
East boundaries. 

An example of such boundary conditions is the domain wall (DW) boundary conditions. For 
the DW boundary conditions the arrows on the boundary of the lattice are going into the lattice 
at the top and bottom of the lattice and are going out of the lattice at the right and left of it. A 
configuration of paths on a 5 x 5 lattice with DW boundary conditions is presented on Fig. [2l 

Notice that the differences 



n{ai)-n{a2), nf, = n(6i) - 71.(62), 



n{C2 



n(ci) 



are the same for all configurations with given fixed boundary conditions. Here n{x) is the total 
number of vertices of type x = ai,bi,Ci in the configuration. 
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Figure 2. A possible configuration of paths on a 5 x 5 square grid for the DW 
boundary conditions 

In particular, the partition function for fixed boundary conditions trivially depends on magnetic 



In this paper we will focus on the 6-vertex model with fixed boundary conditions. 

2.2.2. Another important type of boundary conditions are the periodic boundary conditions. In 
this case the edges at opposite sides of Ln^m are identified so that the configuration of arrows on 
the left and right boundary is the same as well as the configuration of arrows on the top and bottom 
boundary. 

The 6-vertex model with periodic boundary conditions is an example of an "integrable" (solvable) 
model in statistical mechanics and has been studied extensively, see |Baxj . [LW] and references 
therein. In particular, it means that the row-to-row transfer-matrix of the model can be diagonalized 
by the Bethe ansatz. 

2.3. The Height Function. By outer faces we mean unit squares centered at (— ^, m), {N — ^,m), 
with m = — i, i, . . . , M — i and (n, — i), (n, N — i), with n = — i, ^, . . . , N — ^. Each corner outer 
face has two edges in their boundary, other outer faces have three edges in their boundary. 

A height function h is an integer- valued function on the faces -F/v of the grid L^r (including the 
outer faces), which is 

• zero at the southwest corner of the grid, 

• non-decreasing when going up or right, 

• if /i and /2 are neighboring faces, then |/i(/i) — /i(/2)| < 1- 



fields: 








8 



KONSTANTIN PALAMARCHUK AND NICOLAI RESHETIKHIN 

















1 


2 


3 


4 






1 


2 


2 


3 









1 


1 


2 









1 


1 


1 

















Figure 3. Values of the height function for the configuration of paths given on 
figure 2. 



The boundary value of the height function is its restriction to the "outer faces" . Denote the set of 
outer faces by dF]^. Given a function h^^^ on 9-F/v denote ^{{h^^^) the space of all height functions 
with the boundary value h^^^ . 

If we enumerate faces by the coordinates of their centers the height function can be regarded as 
(A^ + 1) X (M + 1) matrix with non-negative entries. 

It is clear that there is a bijection between states of the 6-vertex model with fixed boundary 
conditions and height functions with fixed boundary values. 

Indeed, given a height function consider its "level curves", i.e. paths on the grid L^v, where the 
height function changes its value by 1, see Fig. [3l Clearly, this defines a state for the 6-vertex 
model on L^v with boundary conditions determined by the boundary values of the height function. 

On the other hand, given a state in the 6-vertex model, consider the corresponding configuration 
of paths. It is clear that there is a unique height function whose level curves are these paths and 
which satisfies the condition /i = at the southwest corner. 

It is clear that this correspondence is a bijection. 

There is a natural partial order on the set of height functions with given boundary values. One 
function is bigger then the other if it is entirely above the other. There exist the minimum /imin 
and the maximum /imax height functions such that /imin ^ h < /imax for all height functions h. 

Thus, we can consider the 6-vertex model as a theory of fiuctuating discrete surfaces constrained 
between /imax and /imin- Each surface occurs with probability given by the Boltzmann weights of 
the 6-vertex model. 
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2.4. The inhomogeneous 6-vertex model and volume weights. In the inhomogeneous 6- 
vertex model the Boltzmann weights depend on the edge. Thus, we have 6N^ parameters ai{m, n), 
bi{m,n), and Ci{m,n). 

Let us assume that the inhomogeneity is only in magnetic fields, i.e. weights a,b, and c do not 
change from edge to edge, but the magnetic fields H(e), V{e) do. 

Let {P} be the collection of paths corresponding to a state in the 6-vertex model and {h{f)} be 
the corresponding height function. 

Proposition 2.1. Let us assign weights s(e) to the edges of the lattice and 1 to the outer edges, 
then 

/ P eeP 

where qf = s{ei)s{e2)s{e3)''^s{e4)~^ for inner faces, qf = s(e)^(^) for the outer face f adjacent to 
the edge e, qf = 1 for outer faces not adjacent to any edge of the lattice (comer faces) with £{e) = 1 
for edges at the upper and right sides of the boundary, e{e) = — 1 for edges at the lower and left 
sides of the boundary. 

The proof is an elementary exercise. 

Let s(e) = exp(i7(e) for horizontal edges and s(e) = exp(y(e)) for vertical edges. Then the 
probability of the state with the height function h in such a model is 

D / ^ ^ ^ ^f^^^^ rivertices "^0 (vertex) 
Pq{state) = . 

Zj 

If H{n,m) = H + am and V{n,m) = V + bn, the weights qj are the same for all faces qf = 
exp(a + b) inside the lattice and the probability is given by 

Pistate) = ^'""■'n^^ice. 

z 

where w{vertex) are the Boltzmann weights with constant magnetic fields and 

vom = J2 hif) 

feLN 

is the volume "under" the height function h. 

3. The thermodynamic limit 
3.1. Stabilizing sequence of fixed boundary conditions. 
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3.1.1. Let a = M/N. We place the grid Ljv,m inside of the rectangle D = {{x,y)\ Q < x < 
1,0 < y < a} so that the vertices of the grid are the points with coordinates (7^^, 7\^) where 
n = 1, . . . , AT, m = 1, . . . ,M. 

We recall that a height function is a monotonic integer-valued function on the faces of the grid, 
which satisfies Lipshitz condition (it changes at most by 1 on any two adjacent faces). The height 
function can be regarded as a function on the centers of the faces of the grid, i.e on (n — 1/2, m— 1/2), 
where n = 0, . . . , + 1, m = 0, 1, . . . , M + 1. The points (i, m - i), {N + \,m- \), (n - 5, 5), 
and (n — ^, M + |), where n = 0, ...,iV-l-l, m = 0, ...,M-|-1 correspond to the "outer" faces of 
Ln,m- 

We introduce the normalized height function as a piecewise linear function on the unit square 
with the value 

hT'^{x,y) = ^hN{n,m). 

for IvTT — ^ — ^^"^ j?TT ^ y ^ ^+T- Here hN{n,m) is a height function on Ln^m- Normalized 
height functions are nondecreasing in x and y directions and they satisfy: 

(1) h{x,y)-h{x',y') <x-x' + y-y'. 

if X > x' and y > y'. 

The boundary value of the normalized height function defines a piece-wise constant monotonic 
function on each side of the region which changes by il/AT or do not change between two neigh- 
boring boundary sites. 

Denote the space of such normalized height functions with the boundary value ho by L^^m(^o)- 

There is a natural partial ordering on the set of all normalized height functions with given 
boundary values descending from the partial order on height functions: hi > /12 if hi{x) > h2{x) 
for all x e D. We define the operations 

hiM h2 = minxeD{hi{x),h2{x)), /ii A /i2 = maXxeD{hi{x),h2{x)), 

It is clear that 

hi A h2 > hi, h2 > hi W h2 

It is also clear that in this partial order there is unique minimal and unique maximal height 
functions, which we denote by h^}'^"' and h^"'^, respectively. 

A sequence of stabilizing fixed boundary conditions is a sequence of functions which are 
boundary values of a normalized height function, and which converges as A'' —t- 00 to a function 0o 
on the boundary of D which is non-decreasing along x and y direction, ^(0, 0) = 0, and satisfies 
the condition: 



(2) 



(f){x, y) - (f){x', y') < X - x' + y - y' 
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We will call the function (j)Q the boundary condition for the domain D (or simply the boundary 
condition). Any function on D which satisfies ([2]), which is non-decreasing in x and y directions, 
and coinside with (j)Q at the boundary is called a /leight function on D with boundary condition (pQ. 

Among all possible boundary conditions </>o we will distinguish piecewise linear boundary condi- 
tions with the slope or 1 along the coordinate axes. We will call these boundary values critical. 
It is clear that any boundary conditions can be approximated critical boundary conditions. 

3.2. The thermodynamic limit. If q is fixed and is not equal to 1 in the large volume limit, the 
system will be in a neighborhood of the minimal height function for g < 1 and in the neighborhood 
of the maximal height function for g > 1. One should expect that the partition function and local 
correlation functions will have finite limit. 

When q = exp(-^) for some A, one should expect the existence of the limit shape. We will study 
this limit in Section [71 

3.3. Gibbs measures with fixed slope. 

Definition 3.1. The Gibbs measure of the 6-vertex model on an infinite lattice has the slope {h,v) 
^/ 

h(n + k,m) — h(n,m) 
hm < ^ ■ — ^ ^- — - >= h 

k^oo k 

and 

h(n,m + k) — h(n,m) 
hm < J — - >= V 

fe->oo k 

It is clear from the definition of the height function that the slope should satisfy conditions 
0<h,v<l. 

The slope is simply the average number of horizontal and vertical edges occupied by paths per 
length. 

The important corollary of the result of [Shefj is that, when a gradient Gibbs measure satisfies 
certain convexity conditions, there exists a unique translationally invariant measure. This implies 
that for the 6-vertex model one should expect the uniqueness of such a measure for the generic 
slope. 

Translationally invariant measures can be obtained by taking the thermodynamic limit of the 6- 
vertex model with magnetic fields on a torus. Then the slope is the Legendre conjugate to magnetic 
fields. 

V / \ r ,N{hor). ^. 1 dlogZ 1 
h = hm (fTp) = lim ( ) = lim 1 — 

N,M^oo" ' N,M^oo N ' JV,M^oo 2iVM dH 2 

for a horizontal e and 

1 siogz 1 

V = lim (cJe) = hm ( — — ) = lim , , — -— — 

Af.M->oo^ ' N,M -^oo^ N ' N,M ^oo 2NM dV 2 

for a vertical e. 
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4. The thermodynamic limit of the 6-vertex model for the periodic boundary 

conditions 



The free energy per site in the thermodynamic hmit is 

log{ZN,M) 



f = — hm 



where Zn,m is the partition function with the periodic boundary conditions on the rectangular grid 
Ln,m- It is a function of the Boltzmann weights and magnetic fields. 

Remark 1. Normally, it is expected that the free energy is not identically zero. Physically, this 
means that the "excitations^' have the characteristic length which is much smaller than the charac- 
teristic length of the system. In some cases the free energy is identically zero, then one expects that 
a typical excitation will he comparable with the size of the system. 

For generic H and V the 6-vertex model in the thermodynamic limit has the translationally 
invariant Gibbs measure with the slope {h,v): 

,15/1 19/1 

The parameter 

^_ a^ + b^-c^ 
2ab 

defines many characteristics of the 6-vertex model in the thermodynamic limit. 

4.1. The phase diagram for A > 1. The weights a, b, and c in this region satisfy one of the two 
inequalities, either a>b + coicb>a + c. 

li a > b + c, the Boltzmann weights a, b, and c can be parameterized as 

(4) a = r sinh(A + r]),b = r sinh(A), c = r sinh(r/) 
with A, r/ > 0. 

li a + c < b, the Boltzmann weights can be parameterized as 

(5) a = r sinh(A — r]),b = r sinh(A), c = r sinh(r/) 

with <r] < X. 

For both of these parametrization of weights A = cosh(?7). 

The phase diagram of the model for a > 6 + c (and, therefore, a > 6) is shown on Fig. [5] and for 
b > a + c (and, therefore, a < 6) on Fig. [H 

When magnetic fields {H,V) are in one of the regions Ai,Bi of the phase diagram, the system 
in the thermodynamic limit has the translationally invariant Gibbs measure supported on the 
corresponding frozen configuration. There are four frozen configurations Ai, A2, Bi, and B2, 
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shown on Fig. HI For a finite but large grid the probabihty of any other state is at most of order 
exp(— aA^) for some positive a. 

Local correlation functions are given by the value of the corresponding observable on the frozen 
state 



lim {ae,...(TejN 



where S is the one of the ferromagnetic states Ai,Bi. 

The frozen regions in the {H, y)-plane are described by the set of inequalities. The boundaries 
of these regions can be derived by analyzing the next to the largest eigenvalue of the row-to-row 
transfer matrix. The description is separated into two cases: a > b + c and b > a + c. Notice that 
a ^ b since A > 1. 

• a > b + c, see Fig. [5l 

^i-region: V + H > 0, cosh(2F) < A, 



(e2-f^-6/a)(e'^ -b/a) > [c/af 



„2V 



e^" > b/a, cosh(2ii') > A, 



742-region: 

i^i-region: 
-B2-region: 

b > a + c, see Fig. [6l 

^i-region: 
^2-region: 
Bi-region: 



V + H <0, cosh(2F) < A, 



b/a){e 



-2V 



b/a) > {c/af, 



-2H 



> b/a, cosh(2F) > A, 



(e2^-a/6)(e-2V 
(e-2^^_a/6)(e2^ 



a/b) > {c/bf, 
a/b) > {c/bf. 



^2H 



'2H 



> a/b, 
> a/b. 



i?2-i'egion: 



a/b){e 



[e^^ -b/a){e^^ -b/a) > {c/a)\ 
{e-^^ - 6/a)(e-2^ - b/a) > (c/a)^ 
V- H>0, cosh(2F) < A, 

2V 



> b/a; 



-2H 



> b/a; 



a/b) > {c/by 



^2H 



> a/b. 



V-H <0, 



cosh(2F) < A, 



a/b){e^^ -a/b) > {c/bf 



-2H 



> a/b. 



cosh(2i7) > A, 
cosh(2i?) > A. 



Figure 4. Four frozen configurations of the ferromagnetic phase 
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Figure 5. The phase diagram in the (ff, F)-plane for a = 2, 6 = 1, and c = 0.8 
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Figure 6. The phase diagram in the (ff, y)-plane for a = 1, 6 = 2, and c = 0.8 
The free energy is a hnear function in K and V in the four frozen regions: 

f = -lna-H-V in Ai, 
(6) f = -lnb + H-V in B2, 

f = -Ina + H + V in A2, 
f = -\nb-H + V in Bi. 
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If {H, V) is a region Di or D2 the corresponding translationally invariant Gibbs measure has the 
slope {h, v) given by (l3|). In this phase the system is disordered, which means that local correlation 
functions decay as a power of the distance d{ei, ej) between and ej when d{ei, Cj) — )• 00. 

In the regions Di and D2 the free energy is given by |SY| : 

f{H,V) = min|min|^i-g-(l-2a)y--^ / ln(- - )p{z)dz], 



(8) p{z) = - + — ; dw - — / , 



= 0. 



where p{z) can be found from the integral equation 

z 27ri J(jz — Z2{w) 2iTi Jqz — zi{w 
in which 

^iH = :7T , Z2iw) = + 2A. 

2A — w w 

p{z) satisfies the following normalization condition: 

1 f 
a = / p(z)dz. 

The contour of integration C (in the complex z-plane) is symmetric with respect to the conjugation 
z — )• is dependent on H (see Appendix [B]) and is defined by the condition that the form p(z)dz 
has purely imaginary values on the vectors tangent to C: 

Re{p{z)dz) 

The formula ([7]) for the free energy follows from the Bethe Ansatz diagonalization of the row-to- 
row transfer- matrix. Its derivation is outlined in Appendix [BJ It relies on a number of conjectures 
that are supported by numerical and analytical evidence and in physics are taken for granted. 
However, there is no rigorous proof. 

There are two points where three phases coexist (two frozen and one disordered phase). These 
points are called tricritical. The angle 6 between the boundaries of Di (or D2) at a tricritical point 
is given by 

cos(^) - ^2 + 2 min(a, 6)2 ( A2 - 1) ' 
The existence of such points makes the 6- vertex model (and its degeneration known as the 5-vertex 
model [HWKKj ) remarkably different from dimer models |KOSj where generic singularities in the 
phase diagram are cusps. Physically, the existence of singular points where two curves meet at the 
finite angle manifests the presence of interaction in the 6-vertex model. 

Notice that when A = 1 the phase diagram of the model has a cusp at the point H = V = 
This is the transitional point between the region A > 1 and the region |A| < 1 which is described 
below. 



c2 
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Figure 7. The phase diagram in the {H, y)-plane for a = 1, 6 = 2, and c = 2 



4.2. The phase diagram |A| < 1. In this case, the Boltzmann weights have a convenient param- 
eterization by trigonometric functions. When 1 > A > 1 

a = r sin( A — 7) , 6 = r sin( A) , c = r sin(7) , 

where < 7 < it/2, 7 < A < vr, and A = cos 7. 
When > A > -1 

a = r sin(7 — A) , 6 = r sin( A) , c = r sin(7) , 

where < 7 < it/2, vr — 7 < A < vr, and A = — cos 7. 

The phase diagram of the 6- vertex model with | A| < 1 is shown on Fig. [71 The phases Ai, Bi are 
frozen and identical to the frozen phases for A > 1. The phase D is disordered. For magnetic fields 
{H,V) the Gibbs measure is translationally invariant with the slope {h,v) = {-^-^ff-^-, ~^^^mr~^)- 
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The frozen phases can be described by the foUowing inequahties: 



(9) 



j4i-region: 
^2-region: 



i?2-i'egion: 



-Bi-region: 




-b/a){e 



a/b){e'- 
-a/b){e 



-2V 



,2y 



2V 



b/a) > {c/a)\ 
- b/a) > {c/af 



a/b) > {c/hf 
a/b) > {c/bf 



e^^ > b/a, 




e^^ > a/b, 



e 



2H 



> b/a, 



The free energy function in the frozen regions is stih given by the formulae ([7]). The first 
derivatives of the free energy are continuous at the boundary of frozen phases, The second derivative 
is continuous in the tangent direction at the boundary of frozen phases and is singular in the normal 
direction. 

It is smooth in the disordered region where it is given by ([7]) which, as in case A > 1 involves a 
solution to the integral equation ([8]). The contour of integration in ([8]) is closed for zero magnetic 
fields and, therefore, the equation ([8]) can be solved explicitly by the Fourier transformation [Bax] 

The 6-vertex Gibbs measure with zero magnetic fields converges in the thermodynamic limit to 
the superposition of translationally invariant Gibbs measures with the slope (1/2, 1/2). There are 
two such measures. They correspond to the double degeneracy of the largest eigenvalue of the 
row-to-row transfer-matrix |Bax| . 

There is a very interesting relationship between the 6-vertex model in zero magnetic fields and 
the highest weight representation theory of the corresponding quantum affine algebra. The double 
degeneracy of the Gibbs measure with the slope (1/2, 1/2) corresponds to the fact that there are 
two integrable irreducible representations of sli at level one. Correlation functions in this case can 
be computed using g- vertex operators |JM] . For latest developments see (BJMST] . 

4.3. The phase diagram A < — 1. 

4.3.1. The phase diagram. The Boltzmann weights for these values of A can be conveniently pa- 
rameterized as 

(10) a = r sinh(?7 — A), 6 = r sinh(A), c = r sinh(7/), 

where < A < r/ and A = — cosh rj. 

The Gibbs measure in thermodynamic limit depends on the value of magnetic fields. The phase 
diagram in this case is shown on Fig. [8] for b/a > 1. In the parameterization (llOh this correspond 
to < A < r]/2. When rj/2 < A < r/ the 4-tentacled "amoeba" is tilted in the opposite direction as 
on Fig. m When (H, V) is in one of the Ai, Bi regions in the phase diagram the Gibbs measure is 
supported on the corresponding frozen configuration, see Fig. HI 
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Figure 8. The phase diagram in the {H, y)-plane for a = 1, 6 = 2, and c = 6 



The Ai, Bi regions on the phase diagram are defined by inequahties ([9]). The free energy in these 
regions is hnear and is given by ([7]). 

If {H^ V) is in the region D, the Gibbs measure is the translationahy invariant measure with the 
slope (/i, v) determined by ([3]). The free energy in this case is determined by solutions to the linear 
integral equation ([8]) and is given by the formula 

If {H, V) is in the region A, the Gibbs measure is the superposition of two Gibbs measures 
with the slope (1/2, 1/2). In the limit A — t- — oo these two measures degenerate to two measures 
supported on configurations Ci,C2, respectively, shown on Fig. [H For a finite A the support of 
these measures consists of configurations which differ from Ci and C2 in finitely many places on 
the lattice. 

We notice that any two configurations lying in the support of each of these Gibbs measures can 
be obtained from Ci or C2 via flipping the path at a vertex "up" or "down" as it is shown on Fig. 
[10] finitely many times. It is also clear that it takes infinitely many flips to go from Ci to C2. 

The 6-vertex model in the phase A is disordered and is also noncritical. 
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Here non-criticality means that the local correlation function {af,.af,.) decays as exp(— Q;(i(ej, ej)) 
with some positive a as the distance d{ei,ej) between and Cj increases to infinity. 

The free energy in the j4-region can be explicitly computed by solving the equation ([8]) via the 
Fourier transform |Bax] . 

4.3.2. The antiferromagnetic region. This boundary of the antiferromagnetic region A can be de- 
rived similarly to the boundaries of the ferromagnetic regions Ai and Bi by analyzing next to the 
largest eigenvalue of the row-to-row transfer-matrix. The difference is that for the region A the 
largest eigenvalue will correspond to n = N/2 and to compute it we should use the solution to the 
integral equation ([8]) in the case when the contour of integration is closed. 

This computation was done in |SY| . |LW| . The result is a simple closed curve, which can be 
described parameterically as 

H{s) = E{s), V{s) = E{ri-9o + s), 

where 

= cosh~^ ' 



and 



.dn(fv.|l-^)^' 
\s\ < 2r], 

1 + max(b/a, a/b)e'^ 



max(6/a, a/b) + e'^ 
The parameter u is defined by the equation r]K{v) = nK'^v), where 

fn/2 fn/2 

K{u)= (l-z^sin2(6'))-i/2^0 K'{i^)= (1 - (1 - i/ 
Jo Jo 
The curve is invariant with respect to the reflections {H, V) — t- {—H, —V) and {H, V) 
since the function H satisfies the identities 



i{r]-ip) =E{r] + ip). 





1 










H ■ 














m 


1 r 





Figure 9. Two most probable configurations in the antiferromagnetic phase. 
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up 



down 



Figure 10. The elementary up and down fluctuations in the antiferromagnetic phase. 



This function is also 4r/-periodic: H(4ry + ip) = 

Proposition 4.1. The boundary of the antiferromagnetic region is a real algebraic curve in e^ 
and e^ given by 

, 2 



(11) 



(1 — v cosh^ Vo) cosh^ H + sinh^ Vq — {1 — u) cosh Vq cosh H cosh V 
(1 — z^cosh^ Vo) sinh^ Vq cosh^ y sinh^ H{1 — i/cosh^ H), 



where Vq is the positive value of V on the curve when H = 0. Notice that v depends on the 
Boltzmann weights a,h,c only through rj. 

Proof. The parametric description of the boundary curve implies that 



dn{-{i^-0o + s)\l-u) = ^— 
IT cosh V 



1/^1 N 1 

dn(— s 1 -u) = — 

TT cosh H 

The addition formula for the Jacobi elliptic function dn [ASj 

dnndnu — (1 — m)snncniisnf cnt; 

dn(M + v) = — ^ 2 • 

1 — (1 — mjsn^nsn^f 

can be used to express dnu and dnf for u = K{r] — ^o)/^ and v = Ks/tt in terms of N and V. 
Using elementary identities for the elliptic functions, we obtain 

sn^iicn^u(dn^w — z/)(l — dn^w) = (dnudnw — dn(M + t;)(cn^M + sn^udn^v))^ 

For u = K {i] — 6o) / TT and v = Ks/tt this identity turns into 

(00 O \^000/0 \/ 0\ 

cn ucosh H + sn u — dnu cosh. H coshV \ =cnnsnMcosh y(cosh — 1)(1 — m cosh H). 

Denote dn(-|^(?7 — 6*0)!! — m) = 1/ cosh Vq, then this identity becomes pT]) . □ 
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5. Some asymptotics of the free energy 

5.1. The scaling near the boundary of the D-region. Here we study the asymptotics of the 
free energy of the 6- vertex model as the point (H, V) inside a disordered region approaches its 
boundary on the phase diagram of the model. 

Let us consider the interface between the disordered region D (D2 for A > 1) and the Ai-region, 
see Fig. [7l It is given hy g{H,V) = 0, where 

(12) 9{H, V) = Hb/a + /h^^I/J - 

Let Hq = (Hq, Vq) be a regular point on the interface, i.e. the interface can be parameterized by 
real analytic functions in its neighborhood, then g{HQ, Vq) = 0. 

We denote the normal vector to the interface at Hq by n and the tangent vector by r. A point 
H = (H, V) in the vicinity of Hq can be represented as H{r, s, t) = Hq + r^sn + rtr, where s 
and t are local coordinates in the normal and tangent directions, respectively, and r is a scaling 
factor such that r — )• 0. Let us choose the normal vector n so that it points in the direction of the 
disordered region D, then the point H{r, s, t) belongs to D if s > . 

Theorem 5.1. Let H{r, s,t) be defined as above. The asymptotics of the free energy of the 6-vertex 
model in the limit r — t- is given by 

(13) f{H{r, s, t)) = fiin{H{r, s, t)) + r]{s, t)r^ + ©(r^), 
where fun^H, V) = — ln(a) — H — V and 

(14) r]{s,t) = -K{es + t^f''^ . 

Here the constants k and 9 depend on the Boltzmann weights of the model and on (Hq, Vq) and are 
given by 



and 



'^ = ^djjg{HQ,VQ) 



4 + {dHg{HQ,VQ)f 
2djjg{HQ,VQ) 



where g{H,V) is defined in IT^) . 

Moreover, djjg{HQ,Vo) > and, therefore, 9 > 0. 

We refer the reader to [P] for the details. 



5.2. The scaling in the tentacle. Assume that a > b. The theorem below describes the asymp- 
totic of the free energy function when H — )• +00 and 

(15) i,n(V<.)-^e--<V<lln(Va) + H ^ o.. 
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These values of {H, V) describe points inside the right "tentacle" on the Fig. [5j 
Let us parameterize these values of V as 

1 r2 



where /3 G [—1, 1]. 



Theorem 5.2. When H ^ oo and (3 £ [—1, 1] the asymptotic of the free energy is given by the 
following formula: 

f{H,V) = -hn{ab)-H- 

^e-2^f/3+ iyr^- -/?arccos(/3)) +0(e-4^), 
Zab V vr vr / 

Proof. From the integral equation for p{z) we can derive the large H asymptotics of the density 
function: 

, , / N 1 2Aa , , 

(16) p{z) = - H TT- + \z\ — > oo 

z 

The integration contour is symmetric with respect to complex conjugation z ^ z. The contour is a 
small deformation of the segment of the circle of radius centered at the origin with endpoints 
having arguments ibvra. 

For large H the free energy function is given by 

f{H,V) = mm (-Ina- H - {l-2a)V 

0<Q!<1 V 



27ri Jfj ' a^z — ab' 

As -?/ — )• oo the density is given by (|16|) and, taking into account the asymptotical description of 
the contour of integration we obtain 

f{H,V) = m:m(-hia - H - V + ai2V - Mb/a)) 

a \ 



where 



ab 



e-'^^K(a) 



dw 

2 



c w 
sin(7ra) 



This integral is easy to compute: 

K{a) 

TT 

The minimum occur at 

(17) cos(7rQ) = ^{2V- Hb/a)) e^" . 
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or, at 

cos(7ra) = /3 

The formula (IIGD follows after the substitution of this into the expression for f{H, V). □ 

5.3. The 5- vertex limit. The 5- vertex model can be obtained as the limit of the 6- vertex model 
when A — )• oo. Magnetic fields in this limit behave as follows: 

• a > b + c. In the parameterization ^ after changing variables H = ^ + / , and V = + m 
take the limit 77 — )• cx) keeping A fixed. The weights will converge (up to a common factor) 
to: 

ai : 02 : 61 : 62 : ci : C2 ^ gA+i+m . ^A-Z-m . ^^A _ ^-Xyl-m . g ; l ; l 

• a + c < b. In the parameterization ([5|) after changing variables H = ^ + I, and V = ^ + m 
take the limit 77 — ?■ 00 keeping ^ = A — r/ fixed. The weights will converge (up to a common 
factor) to: 

ai : as : 61 : 62 : ci : C2 ^ (e« - e-«)e'+'" : : e«+'-'" : e«-'+'" : 1 : 1 

The two limits are related by inverting horizontal arrows. From now on we will focus on the 5- vertex 
model obtained by the limit from the 6-vertex one when a > b + c. 

The phase diagram of the 5-vertex model is easier then the one for the 6-vertex model but still 
sufficiently interesting. Perhaps the most interesting feature is that the existence of the tricritical 
point in the phase diagram. 

We will use the parameter 

7 = 

Notice that 7 < 1. 

The frozen regions on the phase diagram of the 5-vertex model, denoted on Fig. [11] as Ai, A2, 
and Bi, can be described by the following inequalities: 

^i-region: m > —I, I < 0, 

e^^" > 1 -7(1 -e~2'), />1; 
(18) A2-region: m < -I, I <0, 

gS^n < 1 _ i(i _ lyl- 

7 



5i-region: (e^' - )(e-"™ - ) > ,^ ' „ , e^' > 



1 — 7 1~7 (1 ~ 7)^ ' 1~7' 
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-1.5 J 



Figure 11. The phase diagram of the 5- vertex model with 7 = 1/4 (/3 = e '^^). 



5.4. The asymptotic of the free energy near the tricritical point in the 5-vertex model. 

The disordered region D near the tricritical point forms a corner 

--/ + 0(/2) < m < -7/ + 0(/2), h^0 + . 

7 

The angle 9 between the boundaries of the disordered region at this point is given by 

27 

cos(6') = — — -. 

l + T 

One can argue that the finiteness of the angle 9 manifests the presence of interaction in the model. 
In comparison, translationary invariant dimer models can only have cusps as such singularities. 

As it follows from results [HWKKj the limit from the 6-vertex model to the 5-vertex model 
commutes with the thermodynamical limit and for the free energy of the 5-vertex model we can 
use the formula 

(19) f{l,m)= lim {F{r]/2 + l,-7]/2 + m) - F{r]/2,-rj/2)), 

where F(H, V) is the free energy of the 6-vertex model. 



Theorem 5.3. Let ^ < k < - and 

I — — ^ 



m = —kl, 
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The asymptotics of the free energy along this ray inside the corner near the tricritical point is given 
by 

(20) /(/, -kl) = ciik, 7)/ + C2ik, + 0{fl^). 
where 

(21) ci(fc,7) = Y^(-(l + A:)(l+7)+4yfc^), 
and 

The proof is computational. The details can be found in jP]. 

The scaling along any ray inside the corner near the tricritical point in the 6- vertex model differ 
from this only by in coefficients. The exponent h^^^ is the same. The details will be given in a 
separate publication. 

5.5. The limit A — ?• — 1~. If m = 1, the region A consists of one point located at the origin. 

2 



When rj — )• 0+ we have A — )• —1— and m — )■ 1 — . Moreover, K' — t- '7r/2, K — )■ and ~ 

1 + ^(1 — m) sin^(ti). Since cosh~"'^(x) ~ :t^J2{af^^\), when x — )• 1+, we have 

\h-a\7] 
a + b 

Since H((/?) is an odd function we obtain the following asymptotic of H |LW| : 

2 

H((/?) ~ 4e sin( — (/?). 

2r/ 

In this limit the antiferromagnetic region degenerates into the origin H = and V = exponentially 
fast. We note that the point H = V = is special for |A| < 1. 



5.6. The convexity. The following identity holds in the region D [NKj : 



(23) fH,Hfv,v - fuy = \ — ) 



Here g = The constant Dq does not vanish in the L>-region including its boundary. It is 

determined by the solution to the integral equation for the density p{z){see Appendix C). 
Directly from the definition of the free energy we have 

<(n(L)-n(«y> 

N,M^oo NM 

where n{l) and n(r) are the number of arrows pointing to the left and the number of arrows pointing 
to the right, respectively. 

Therefore, the matrix didjf of second derivatives with respect to H and V is positive definite. 
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As it follows from the asymptotical behavior of the free energy near the boundary of the D-phase, 
despite the fact that the Hessian is nonzero and finite at the boundary of the interface, the second 
derivative of the free energy in the transversal direction at a generic point of the interface develops 
a singularity. 

6. The Legendre Transform of the Free Energy 

The Legendre transform of the free energy 

sup(xH + yV + f{H, V)) 

as a function of (x, y) is defined for — 1 < y < 1. 

The variables x and y are known as polarizations and are related to the slope of the Gibbs 
measure as x = 2/i — 1 and y = 2v — \. We will write the Legendre transform of the free energy as 
a function of (/i, v) 

(24) G[h, v) = sup ({2h -l)H + {2v-l)V + f{H,V)). 

a{h, v) is defined on < h,v < 1. 

For the periodic boundary conditions the surface tension function has the following symmetries: 

a(x, y) = a{y, x) = a{-x, ~y) = a{-y, -x). 

The last two equalities follow from the fact that if all arrows arc reversed, a is the same, but the 
signs of X and y arc changed. It follows that afi{h,v) = ay{v,h) and ay{h,v) = afi{v,h). 

The function f{H, V) is linear in the domains that correspond to conic and corner singularities 
of a. Outside of these domains (in the disordered domain D) we have 

(25) VaoV/ = idD, V/oVc7 = idv/(i5). 

Here we gradient of a function as a mapping ^ R^. 

When the 6- vertex model is formulated in terms of the height function, the Legendre transform 
of the free energy can be regarded as a surface tension. The surface in this terminology is the graph 
of of the height function. 

6.1. Now let us describe some analytical properties of the function a{h,v) is obtained as the 
Legendre transform of the free energy. The Legendre transform maps the regions where the free 
energy is linear with the slope (ibl,ibl) to the corners of the unit square V = {{h,v)\ < h < 
1,0 <v < 1}. For example, the region Ai is mapped to the corner h = 1 and v = 1 and the region 
Bi is mapped to the corner h = 1 and v = 0. The Legendre transform maps the tentacles of the 
disordered region to the regions adjacent to the boundary of the unit square. For example, the 
tentacle between Ai and Bi frozen regions is mapped into a neighborhood of /i = 1 boundary of 
V, i.e. 1 and < v < 1. 
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Applying the Legendre transform to asymptotics of the free energy in the tentacle between Ai 
and Bi frozen regions we get 

H{h, v) = -^ln (^—1^^) , V{h, v) = ^ Hb/a) + ^{l-h) cot(^(l - v)), 

and 

(26) a{h,v) = {I- h) In ( '''f . /"^ .. ) - (1 - /i) + t. ln(6/a) - ln(6), 

\ sin(7r(l — V)) J 

Here /i — )• 1— and < f < 1. From ([2^]) we see that cr{l,v) = vln{b/a) — ln(6), i.e. a is linear on 
the boundary h = 1 of D. Therefore, its asymptotics near the boundary /i = 1 is given by 

a{h, v) = v ln(6/a) - In(fe) + {1 - h) ln(l - h) + 0(1 - h), 

as /i — 7- 1— and < v < 1. We note that this expansion is valid when (1 — h)/ sin(7r(l — v)) <C 1. 

Similarly, considering other tentacles of the region D, we conclude that the surface tension 
function is linear on the boundary of P. 

6.2. Next let us find the asymptotics of a at the corners of D in the case when all points of 
the interfaces between frozen and disordered regions are regular, i.e. when A < 1. We use the 
asymptotics of the free energy near the interface between Ai and D regions (I13p . 

First let us fix the point {Hq,Vo) on the interface and the scaling factor r in (jl3p . Then from 
the Legendre transform we get 

and 

3 K(^S + t2)l/2 

It follows that 

1-h edug + 4rt 



l-v 2{-9 + rtdHg)' 
In the vicinity of the boundary r — t- and, hence, 

1^ ^ _dHg ^ 1 - b/a e-^^° 
^ ' l-v 2 l-b/ae-^^o 

as /i, f — 7- 1. Thus, under the Legendre transform, the slope of the line which approaches the corner 
h = V = 1 depends on the boundary point on the interface between the frozen and disordered 
regions. 

It follows that the first terms of the asymptotics of a at the corner h = v = 1 are given by 
a{h,v) = -Ina - 2(1 - h)Hoih,v) - 2(1 - v)Voih,v), 
where Ho{h,v) and Vo{h,v) can be found from ([271) ^^id g{Ho,Vo) = 0. 
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When |A| < 1 the function a is strictly convex and smooth for sdl < h,v < 1. It develops 
conical singularities near the boundary. 

When A < —1, in addition to the singularities on the boundary, a has a conical singularity at 
the point (1/2, 1/2). It corresponds to the "central flat part" of the free energy /, see Fig. [HI 

When A > 1 the function a has corner singularities along the boundary as in the other cases. 
In addition to this, it has a corner singularity along the diagonal v = h if a > b and v = 1 — h if 
a < b. We refer the reader to [BSj for further details on singularities of a in the case when A > 1. 

7. The Thermodynamic Limit and the Variational Principle for Fixed Boundary 

Conditions 

7.1. The Variational Principle. 

7.1.1. Let a be the surface tension function of the 6- vertex model with the periodic boundary 
conditions defined in (|24p . We consider the functional 

(28) I[ip] = I a{Vip)(fx + X I ip(fx, 

Jd Jd 

Let h(x,y) be a minimizer of this functional on the space L{D,ipQ) of functions nondecreasing 
in X and y directions and satisfying the condition 

^{x, y) - ^{x, y') < X - x' + y - y' 

and the boundary condition 

'^\dD = ^0- 

Notice that height functions are Lipschitz with \ip{x,y) — ip{x',y')\ < 2max{\x — x'\, \y — y'\). 

Proposition 7.1. The functional I[ip] has unique minimizer. 

Indeed, since a is convex, the minimizer is unique when it exists. The existence of the minimizer 
follows from compactness of the space L{D,(Pq) in the sup norm. The arguments are completely 
parallel to those in [CKPj. 

7.1.2. If the vector V/i(a;,y) is not a singular point of a, the minimizer h satisfies the Euler- 
Lagrange equation in a neighborhood of (x, y) 

(29) div(Vcj o V/i) = A. 
We can also rewrite this equation in the form 

(30) Va{\/h{x,y)) = ^ {x,y) + {-gy{x,y), g^ix^y)), 

where g is an unknown function such that gxy{x,y) = gyxix,y)- It is determined by the boundary 
conditions for h. 
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Applying ([25]) . we see that 

(31) yh{x,y) = Vf(^^x - gy{x,y),^y + g^{x,y)y 

From the definition of the slope, see ^ we have < /h < 1 and < /y < 1. Thus, if the 
minimizer h is differentiable at {x,y), it satisfies the constrains < hx < 1 and < hy < 1. 

7.2. Large deviations. The following statement is a minor variation of the theorem 4.3 from 

Theorem 7.2. Let N — )• oo, A be finite and q = exp(-^) then the sequence of random normalized 
height functions h^ converges in probability to the minimizer of 128\). The rate of convergence is 
exponential of N"^ . 

The minimizer of the variational problem psp is called the limit shape of the height function. 

This theorem is the manifestation of the general philosophy of the large deviations principle. 
The probability of a having state with the height function h has the has the following asymptotic 
N ^ oo : 

Prob{h) ~ exp(^AiV2 j hd^x + j (T{^h)d^i 

Here a is the surface tension function for the periodic boundary conditions. Clearly this probability 
has maximum at the limit shape. States with the height function, which macroscopically differ from 
the limit shape, should be expected to be exponentially improbable. The theorem states that this 
is exactly what is taking place. 

8. The limit |A| — ^ oo 

8.1. Minimal and maximal height functions. The space of normalized height functions on 
L]\j^M has the partial order described in section 2. Denote the minimal and maximal normalized 
height functiions with respect to this partial order /i™" and /i""*^, respectively. 
For two functions hi, /12 on D define the distance 

(32) dist{hi, h2) = supxeD\hi{x) - h2{x)\ 

Similarly define the the distance between two functions on the boundary of D. 

Let h^^ be a sequence of functions on the boundary of D converging to ipo, and such that h^^^ 
is a boundary of a normalized height function of Ln,m- Denote /i^^ and hl^Jx the minimum and 
maximum height functions from L7v,m(^o^^)- The following is clear: 

Proposition 8.1. Let ipmin o-nd (pmax be minimal and maximal height functions with the boundary 
condition ipQ. Then, <pmin = hmjv_>.oo /'-jv*" ^^'^ Vmax = hmjv_>.oo ^tv'*^ with respect to the distance 

m 
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The functions (pmin and (p. 



max 



minimize the functionals 




Jd 



in the sapec H{ipo). 

Let us assume that the boundary conditions are critical, that is y?o is piece-wise hnear, non- 
decreasing in X and y direction with the slope or 1. In this case the minimum and maximum 
height functions are piecewise linear functions, such that each linear part has the slope either or 
±1 along coordinate axes. 

We will say a point {x, y) is regular if at this point the function ipmin is differentiable. For 
critical boundary conditions regular points form regions with piece-wise linear boundary where the 
function (pmin has a constant slope. We will call them linear domains. 

Points where the gradient of ipmin is discontinuous will be called singular points. For critical 
boundary conditions singular points are the points where several linear domains meet. The ualency 
of a singular point is the number of linear domains which meet at this point. For generic critical 
boundary conditions the valency of each critical point is at most three. 

The list of all possible phases at a tricritical point of the 6-vertex model with generic critical 
boundary conditions is given on Fig. [I2j 

8.2. The asymptotic of the minimizer when |A| — >• oo. Here we study the asymptotic of 
minimizer of I\[<f] as A — )• ±00. It is more convenient to divide I\ by |A|, so we are looking for the 
asymptotics of the minimizer hx of 



Let us focus on the limit A — )• -|-oo. The limit A — )• —00 can be treated similarly. 

When A — )• -|-oo, the minimizer hx approaches the minimal height function (pmin described above. 
Let us look for the asymptotical formula for the minimizer in a small neighborhood of a point (xq, 2/0) 
of the form 

hx{x,y) = (pmin{x,y) + J Hx{X{x - xo), A(y - yo)), 
where Hx = H + o(l) as A — )• -l-oo. 

Let z = (2:1,22) = {r,s) = (A(x - xo),A(y - yo)) and x = (xi,X2) = {x,y). Because (pmin is 
linear, its second derivatives vanish and we can rewrite the Euler-Lagrange equation ()29p as 



Notice that for critical boundary conditions the function V(/?mm is piece-wise constant. 

Now assume that (xo,yo) is a singular point, i.e. a point where two or more linear domains of 
i^Pmin meet. Recall that for generic critical boundary conditions only two or three linear domains 
can meet at a point. Taking the limit A — )• 00 in (|33p we obtain. 




(33) 
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Figure 12. Possible phases at a tricritical point of the 6- vertex model. 

Proposition 8.2. Let h\ he the minimizer of and (xo,yo) be a singular point, then for each 
(s, t) G M? there exits 

(34) H{r, s) = lim A (hxixo + ^,yo + ^)- (^™"(xo, Vo)) 
This function is the solution to 

(35) dw{Vao{VH)) = 1, 
with the boundary conditions 

(36) H{r, s) hr + ks 

as (r, s) — )• oo along a ray in the i-th linear domain adjacent to (a;o,yo)- Here {ki,li) = Vipmin in 
the i-th linear domain and since the boundary conditions are critical, ki,li = 0, 1. 

Conjecture 8.1. The function H is once differentiable. There is a smooth curve separating 
into regions where H is linear and regions where H is smooth with positive definite matrix of second 
derivatives. 
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Remark 2. For non- generic boundary conditions in the thermodynamical limit more then three 
linear domains can met at one point. In this case one should expect that the conjecture still holds 
with more then three linear domains meeting at a singular point. 

8.3. The asymptotic near double degenerate singular points. 

8.3.1. A height function defines the surface z = h{x, y) in M^. Regions where h is linear are planes. 
Here we will describe the solution to the equation (j35p with the boundary conditions (j36p in the 
case when there are only two asymptotic planes meeting at(xo,yo) • 

Let {ki, and {hjh) be the directions of the steepest assent of these planes. The numbers ki 
and li are either or 1 since we assume critical boundary conditions and therefore (fmin is piece- wise 
linear with slope (0, 0), (0, 1), (1, 0) or (1, 1). 

The function H{r,s), defined in (j34l) . has the asymptotic conditions (|36p determined by these 
planes. The function H{r, s) is also invariant with respect to translations in m = (/c — Z)-'--direction. 

Thus, we are looking for a function n such that 



with some constant C. 

8.3.2. For the domain wall boundary conditions minimal and maximal height functions are shown 
on Fig. [131 In this case the boundary between two linear domains is a line in the direction (1, — 1) 
for the minimal height function and is the line in the direction (1, 1) for the maximal height function. 

Thus, for these boundary conditions the function H describing the asymptotic of the minimizer is 
constant in the direction (1, —1) when A — )• +00 and it is constant in the (1, l)-direction if A — t- —00. 

Let H{r, s) = k{s + r) be the solution to (|37p which is invariant with respect to translations in 
the (1,-1) direction. The symmetries of the Legendre transform a of the free energy imply that 
ah{h,v) = av{v,h). Using the equation (1381) and this symmetry of a, we obtain 



H{r, s) = K {{ki - h)r + (fca - h)'?) , 

which satisfies the differential equation (j35p with the asymptotic conditions ([36]). 
Let us introduce 




Then the differential equation for H becomes the first order ODE for the function k' 
(37) K"{t)S{K'{t)) = I. 

Integrating it, we obtain the equation defining the function g{t) = n'{t) implicitly 



(38) 




— {g{t),g{t))=t/2 + C/2. 
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Taking into account (j25p we obtain 

g{t) = dHf{t/2 + C/2, t/2 + C/2) = dvf{t/2 + C/2, t/2 + C/2) 

and, hence, 

Kit) = f{t/2 + C/2,t/2 + C/2). 

In the case of two asymptotic planes we have an additional symmetry of H{r, s) with respect to 
the intersection line of these planes, i.e. K{t) = K{—t). The free energy / also has the symmetry 
f{x,y) = f{—y, —x). Therefore C = and we proved the following statement. 

Theorem 8.3. For DW boundary conditions the minimizer h\ when A — )• +oo has the following 
asymptotic when x = + r / \, y = —xq + s/X: 

(39) hx{x,y) = ^min + ^ f ,^^\ +o[j^ 
If x = xq + r / X and y = xq + s / X and A — — oo the asymptotic of the minimizer h\ is given by 

1 I T — S T — 5\ /I 

(40) hx{x,y) =if.ratn + ^f , 2~ ) ^ "K)? 

8.3.3. When c> a + b (which is equivalent to A < 0) the height function hx develops extra linear 
domains known as facets. These are the regions of the antiferroelectric phase discussed in section 
14.3.21 In dimer models this phenomenon is studied in |KOS| . 

The facets also appear in the function H describing the asymptotic of the minimizer as A — )• oo. 

For the function H describing the asymptotic near a double singular point (xcyo) where two 
linear domains meet along the diagonal x + y = the facet is a strip 



in coordinates x = xq + f, y = —xq + j- Its width R is given by the formula 



R = V2 -(^ 



h, 



'mm 



Figure 13. The minimum and maximum height functions for the DW boundary conditions. 



34 



KONSTANTIN PALAMARCHUK AND NICOLAI RESHETIKHIN 



In the hmit A— t-— 1— orr/— )• 0+ the asymptotic of H gives the foUowing asymptotical value of 

R: 



_2 



(41) i? = 4^/2e 2. sin(^^-^^j(l + 0(7?)). 

The free energy / is linear in the antiferromagnetic region. It is growing with exponent 3/2 in 
the normal direction to the boundary outside of the antiferromagnetic region. This agrees with the 
Pokrovsky-Talapov law which states that / should be growing with exponent 3/2 in the normal 
direction to the boundary of the facet |PTj . In particular, h grows with the exponent 3/2 in the 
(1, l)-direction near the boundary of the c-droplet far enough from the boundary of the square. For 
s — )• + we have 

hxixQ + ^, -xo) = /lo + ^(s - ^)^/^ 
where /iq is the value of the height function at the boundary of the facet and «: is a constant which 
can be computed explicitly. 

9. Conclusion 

9.1. Correlation functions in the bulk. As we have seen in the previous section at the macro- 
scopical distances in the thermodynamical limit the height function is deterministic and is the 
minimizer for the variational problem (|28p. 

But the height function at smaller distances remain random. Their fluctuations are described by 
the asymptotical behavior of correlation functions in the thermodynamical limit. These asymtotics 
have been studied extensively in dimer models which describe the A = case of the 6- vertex model. 

The thermodynamical limit of correlations functions in the bulk describes translationary invariant 
Gibbs measure with given polarization. These asymptotics of correlation functions have been 
studied a lot using various methods which are essentially based on representation theory of affine 
quantum groups. 

The exact computation of the asymptotic of correlation functions for local observable separated 
by large distances on the lattice remain one of the main problems. On the other hand this is one of 
the most important physically relevant information about correlation functions . Some information 
about this asymptotic of correlation functions can be obtained using the arguments of finite-size 
scaling and the assumption of conformal invariance of the leading terms of the asymptotic. 

Let us consider the 6- vertex configurations of paths which may end at some edges. The weight 
of such configurations is given the same product of weights as before. Define local observables Te 
and T* as 

1, if a path going up starts at e; 
0, otherwise 



r:iS) = 



1, if a path going down starts at e; 
0, otherwise 
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The value of a product of such observable when each of the factors correspond to a different edge 
is the product of values of observables. 

The following formulae were obtained in |BIRj for H = and when all edges are vertical on the 
same row using the finite-size scaling and the assumption of conformal invariance: 

, , , 9 ^ Bcos(2kFd) 

(42) <ae,ae,>c^<ae>' +-^ + ^ '- + ... 

(43) <r;^re, >~ -^ + ... 

Here j is the distance between ei and 62, a = 27r^p(^)^. The terms denoted by . . . are with higher 
powers of d~'^ and d~°' where d is the distance between ei and 62- These higher order terms may 
also be oscillating. Because H = the integration contour C in the integral equation fro p{z) is a 
segment of the real line in the additive parameterization. The constant kp is the Fermi-momentum 
and it is also can be expressed in term of p{z). It is also equal to the vertical electrical polarization. 



The finite size computations were extended to the case V 7^ in jNKj . They argued that the 
complete spectrum of effective c = 1 conformal field theory is given by 



1 TnP 

A± = — ( h n^g lb 2nm), n, m G Z 

4 g 

where g is defined by the Hessian of the free energy and is related to a as a = 

Combining these results we conjecture that for generic electric fields the asymptotical behavior 

of correlation functions is still given by formulae (j42p and (|43p with a = ^ 

One can show easily that when A = we have g = 1/4 and the asymptotics (|42|) (^51 agrees with 

the results from IK I on dimer models. 



9.2. Open problems. Here we will list some open question about the 6-vertex model and other 
related models. 

• Find the classification of generic singularities of limit shapes. For A > 1 one should expect 
the presence of corners as generic singularities of limit shapes. 

• Find the scale at which fluctuations near singularities of limits shapes are described by some 
random process and describe such processes. For example in dimer models such fluctuations 
near the boundary of the limit shape are described by Airy process, and near a generic cusp 
are described by the Pearcey process. 

One can argue that the same processes should describe correlation functions near similar 
singularities of the limit shape in the 6-vertex model but this is still a conjecture. 

The scaling near the corner singularity seems particularly interesting problem since such 
singularities do not appear in dimer models. 
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• Understand the role of integrability of the 6-vertex model in the formation of the limit 
shape. Here by integrability we mean that the model can be solved by the Bethe ansatz, 
and that the weights satisfy the Yang-Baxter equation (and therefore transfer-matrices form 
a commuting family). 

• The 6-vertex model is closely related to the representation theory of quantized universal 
enveloping algebra of sl2- It would be extremely interesting to see which aspect of the 
representation theory of this algebra naturally appear in the limit shape phenomenon and 
in the scaling of correlation functions near singularities. 

• The 6-vertex model has natural generalizations related to other simple Lie algebras. Con- 
figurations in these models can be described in terms of r height functions where r is the 
rank of the Lie algebra. The limit shape in this case is a surface in M'^"'"^. It would be 
extremely interesting to investigate such systems trying to answer all questions mentioned 
above. 
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Appendix A. The Partition Function of the Inhomogeneous 6- Vertex Model 

Here we compute the partition function of the inhomogeneous 6- vertex model with the periodic 
boundary conditions applying the Bethe Ansatz method. We briefly describe the results and refer 
the reader to [KBIj for the details about the method. 
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Bethe Ansatz method works for the inhomogeneous 6-vertex model with the fohwing inhomo- 
geneities. One should parameters Ai, A2,.-) to each row of the grid Lat^m and parameters fii, 
H2t-, fJ-N to each column of the grid as shown on Fig. [T^ Thus, a pair of spectral parameters 
(Aj,/ij) is associated to the vertex of the grid. Boltzmann weights assigned to this point are: 



ai{i,j) = a{Xi - ij,j)e"~^^, ai{i,j) = a{Xi - ^j)e 



H+v 



-H+v 



hihj) = b{^i - fJ'j)e ^ , bi{i,j) = b{Xi - Hj)e~ 

ci{i,j) = c{Xi - fij), ci(i, j) = c{Xi - fij) 

where functions o(A), &(A), c(A) describe the parameterization of Boltzmann weights of the homo- 
geneous model. 

The partition function of the model with periodical boundary conditions can be written as 

N 

ZN = Ttl[T{Xk), 

k=l 

where T is the 2^ x 2^ row-to-row transfer matrix of the 6-vertex model |Bax| . The raw-to-raw 
transfer- matrix is the trace of the " quantum monodromy matrix" : 

r(Afc) = tror(AA:). 

Here r is the product 

r(Afe) = i?o,i(Afc - /ii) . . . Ro^Ni^k - IJ-n), 
which acts in Cg) C^®^. The matrix RQj(Xk — fJ-j) acts as the matrix 

fai{Xk-fij) \ 

b2iXk-Hj) ci(Afc-//j) 

C2(Afc-/ij) hi{Xk-lJLj) 

V a2{Xk-Hj)J 

in the basis ei (gi ei, ei (8) 62, 62 ^ ei, 62 (X" 62, of the product of the 0-th and j-th factors in the tensor 
product. It acts trivially in other factors. 



^1 [i-l ^3 



Figure 14. The 4-by-4 square grid with corresponding spectral parameters. 
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Let us denote the number of arrows pointing down in a row of the grid by n. Denote by (C®^)„ 
the corresponding subspace in the space of all possible states on the raw of length N of vertical 
edges. The action of the transfer- matrix preserve these subspaces. 

The Bethe Ansatz method gives the following result for the eigenvalues of T: 

(44) A(A,ui,..,u„,i/,y) = KL{Kui,..,Un,H,V) + KR{\,ui,..,Un,H,V), 
where 

n , . N 

(45) A.(A, u„..,n^,H,V) = e-^^+(^-2«)V J] JJ 6(A - ^,), 

1=1 ^ j=i 

and Ui satisfy the Bethe equations 

/^gs ^2NH TT ajUk - Ul)b{ui - Uk) ^ jr b{ui - ^j) 

k=il^i ~ Uk)b{uk - ui) a{ui - Hj) ' 

It follows that the partition function of the model is given by 

N 



^ Aa(Afc,lil,li2, ..,Un,H,V). 



a k=l 

where Aq are the eigenvalues of the row-to-row transfer matrix T and the sum is taken over all 
eigenvalues with their multiplicities. 

Thus, in the homogeneous case when all are the same the asymptotic of the partition function 

as M — )• cxD is 

where d is the multiplicity of the largest eigenvalue h.max{N) of the row-to-row transfer-matrix on 
N sites . 

Therefore, if we can compute the asymptotic of the largest eigenvalue of the transfer-matrix as 
— )• oo we can find the asymptotic of the partition function in the thermodynamical limit. 
Strictly speaking, this logic gives the asymptotic of the partition function in the limit M S> A S> 
1. Under mild assumptions one can argue that the leading term of the asymptotic of free energy is 
uniform. However, in some cases a resonant phenomenon may occur which will make the leading 
term of the asymptotic of Z^^m depend on the ratio N/M |KW| . 
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Appendix B. The Free Energy of the Homogeneous 6- Vertex Model 

In the homogeneous 6-vertex model the Boltzmann weights are the same for all vertices on the 
grid Ljv- When | A < 1 the functions a, b, c are: 

a = p s'm{u + r]) , b = psm{u), c = psin[r]). 

Other values of A can b obtained by analytical continuation. 

The formula for eigenvalues of the row-to-row transfer matrix given in the previous section 
becomes 

A(")('iV) = a^Q^H-\-{N-2n)V TT sin(-Uj ~ u + rj) ^ ^-NH+{N-2n)V TT sin(ttj - u - rj) 

siniuj — u) sm(un — u) ' 

j=i J j=i J 

where {ui, ..jUn) are solutions of the Bethe ansatz equations 



N 

sin(nj) 



sin(nj -|- rj) 
Introducing new variables 



2NH A sin(txfc - Uj + rj) 
^J-^ sm{uk - Uj - rj) 



sinf^u) 

-, j = l,..,n, w- 



sm(uj + rj) sm{u + rj) 

we get 

I&^HN) = (j-'VgiVi?+(Ar-2n)y TT 1 - 2Aw + WZj ^ ^N^-NH+{N-2n)V TT ^ 2Azi + WZj 

J. J. — in J-J- 



Zj — W -*■ W — Zj 

where Zi satisfy the Bethe equations 

(47) zf = (-l)«-ie2^^ ft l-2A^j+^i^fc 
^ ' 3 y ) ni_2Azk + ZjZk 

Wc want to find the asymptotic of the largest eigenvalue of the transfer-matrix when AT — >■ oo. 
We shall do it in two steps. First we will find the asymptotic of the largest eigenvalue among 
A^*") (N) in the limit 

n 

(48) N — > oo, n — >■ oo, ~ ^ ~ const, 
After this we will find for which < a < 1 this eigenvalue is the largest. 

Appendix C. The 6- Vertex Model Bethe Equations 



Let us denote 

(49) 



^ie{z,w) ^ 1 - 2Az + zw 
1 — 2Aw + zw 
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Then, taking the logarithm of both sides of each Bethe equation ()47p and dividing by A^, we get 

n — 1 i " 

(50) ln{zj) = ':^j^iTT + 2H + ^Y,^(^j^^k) + 2TTikj/N, j = l,..,n, 

k=l 

where kj are integers. 

All theoretical and numerical results on the 6-vertex model suggest that the following is true. 

Conjecture C.l. In the limit N oo with — l/S fixed the roots of the Bethe equations 
corresponding to the maximal eigenvalue of the transfer matrix in the subspace (C®^)„ concentrate 
on a smooth curve C in the complex z-plane. The curve C is symmetric with respect to the complex 
conjugation z ^ z. 



sponds to kj — — 



Conjecture C.2. The maximal eigenvalue of the transfer-matrix in the subspace (C® )„ corre- 
— J- 



These conjectures can be proven in the free-fermion case A = 0. 

The conjecture (IC.ip implies that as A'^ — t- oo one should expect finite limit density of zj on the 
contour C. li Zj ^ z £ C the limit p(z) = liuiiy^oo Wf — — r is the value of the limit density 
function of solutions to the Bethe equations on C. Let us denote the endpoints of the contour C 
by ^ and ^, ^ = limiv_^oo zi and ^ = limAr^oo Zn- 

Taking the limit in the equation (I50p we obtain the integral equation for p{z): 

(51) ln{z) = 2H — TTia -\ / 9{z,w)p{w)dw + [ p{w)dw. 

27r Jc 

Here z € C. Differentiating this equation with respect to z, we obtain another, more convenient 
form of this equation: 

(52) p{z) = i + -L / K{z,w)p{w)dw, 

z 2711 Jc 

where 

(53) K{z, w) = -zdj{z, w) = ""^ + ^ ^7 ^ . 

1 — 2Az + zvu I — 2Auj + zw 

The kernel K is singular on two curves 



(54) 



Ci = {z = zi{w) = weC}, 



C2 = {z = Z2{w) = -- +2A\ weC}. 
w 



and it can be written as 

(55) K{z,w) = + 



Z — Zi{w) Z — Z2{w) 

Notice that functions Z2{w) = z^^{w), i.e. zi and Z2 are inverse to each other. 
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If A = 0, coutours Ci and C2 coincide and the kernel is zero. 

Conjecture C.3. //A 7^ 0, there are two cases: 

• Conours C,Ci, C2 do not intersect each other; 

• Contours C,Ci, C2, all intersect at two points w± = e^*'' which are the solutions of up' — 
2/\w+l = and are conjugate to each other (these two points may coincide in the degenerate 
case). This is possible only if a = 1. 

The integral equation on p{z) can be rewritten as 
/^^\ / N 1 1 /" Piw) , 1 f p(w) 



z 27rz J(jz — Z2{w) 2-Ki J(jz — zi{uj) 

Notice that the contour C can be deformed as long as it does not intersect the curves Ci and C2. 

The contour C is defined by the condition that the form p{z)dz has purely imaginary values on 
the vectors tangent to C 



= 0, 

zee 



(57) Re{p{z)dz 
and by the normalization condition on p{z) 

(58) ^ I p{z)dz = a. 
The free energy per site is defined as 

The formula for eigenvaues A*^"^ {N) has two terms. Each of them grow exponentially for large N . 
Generically, one of the terms will dominate the other. Taking into account conjectures about roots 
of Bethe equations in the thermodynamicallimit we obtain the follwoing integral representation for 
the free energy [LWj : 

r r I [ w'^ - 2Aw + 1 ^ 
f = mini mini El — H — (1 — 2a)V / ln(w )p(z)dz>, 

I a I 2'Ki Jq W — Z ) 

min\E2 + H -{l-2a)V - — [ ln{2A - w + ~ '^^'^ ^ ^ )p{z)dz}\, 
a L 2m Jc w — z ) ) 



or 

„2 



f{H,V) = min|min|gi-ij--(l-2a)y--^ / ln(- - / , )p{z)dz}, 

min|s2+i?-(l-2a)l/--iT /" ln(^^^ + — ^^)p(z)dz| |, 
a L I'm J (J ab ab — a^z i ) 

where p{z) is determined bythe integral equation ([52]) with conditions ([57]) and ([58]) . 
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Appendix D. The Free Energy in the Free-Fermion Case 



The free-fermion curve in the space of parameters of the 6-vertex model is A = 0. Since in this 
case K{z,w) = the integral equation of p immediately implies p{z) = 1/z. The contour C in this 
case is an arc 



(59) 



where 6 £ [—TTa,-Ka]. 

Let us assume that a > b. The formula (I59p for the free energy with fixed polarization a gives 



(60) 



fa = Ei-H-il- 2a)V - ^ n In 



a' 



+ 62 



a ab — a^e^^e*^ 



Minimizing this expression in a we get the equation 



(61) 



cos(7ra) 



„2(g2//+2y _g-2^^-2l/) + 52(g- 



2_fi"+2y _ g2^i"-2y^ 



This defines the critical value > a > 1 when the absolute value of the r.h.s is not greater then 1. 
Otherwise a = or a = 1. 



Proposition D.l. The free energy can be written as the following double integral: 



f = mino<a<ifa 



2n 



2n 



Proof. The proof is computational. For the double integral in question we have: 



where 



J de J dipln(^ai + 6ie^^ + (62 - a2e^^)e^^) 
2tt / dein{ai + bie''^) + 2-K / ^6* ln(62 + 026 



if \a2 - 62I > fli + bi, 

TT if |ai - 6i| > 02 + 62, 
_ _ a2 + 62 _ ^2 _ ^2 



cos 



otherwise. 



2(ai6i + 0262) 

The last equality holds if |a2 — 62I < ai + 61 and |ai — 6i| < 02 + 62. 
We can rewrite / as 



(62) 



8^2/ = 8vr2lnf^i±^l±i^l^ 



) +2tt [ dein 


62 + 026'^ 




ai - 6ie*^ 
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We note that the argument of the logarithm under the integral is greater than 1 for any 9 in the 
domain of integration. 

Taking into account the parameterization of Boltzmann weights by electric fields H and V it is 
easy to see that (1620 together with the equation for 9* coincide with the formula for the free energy 
derived from the Bethe Ansatz. 

□ 

The double integral formula can be obtained directly from the Pfaffian solution of the dimer 
model which maps to the 6-vertex model at A = 0. 
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